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The coefficient of normal restitution of colliding viscoelastic spheres is computed as a function 
of the material properties and the impact velocity. From simple arguments it becomes clear that 
in a collision of purely repulsively interacting particles, the particles loose contact slightly before 
the distance of the centers of the spheres reaches the sum of the radii, that is, the particles recover 
their shape only after they lose contact with their collision partner. This effect was neglected in 
earlier calculations which leads erroneously to attractive forces and, thus, to an underestimation of 
the coefficient of restitution. As a result we find a novel dependence of the coefficient of restitution 
on the impact rate. 

PACS numbers: 45.70.-n,45.70.Qj,47.20.-k 



I. INTRODUCTION 

The dynamics of a granular system is governed by the 
particle interaction law, that is, by the forces the parti- 
cles in contact exert on one another. In general, these 
forces may be complicated functions of the time depen- 
dent mutual deformation and relative velocities in normal 
and tangential direction. In the simplest case the parti- 
cles are modeled as spheres interacting via normal and 
tangential forces. 

Given particles of radii R\/2 and masses m 1 / 2 at po- 
sitions fi(t) and ^(i) traveling at velocities vi(t) and 
V2(t). The particle deformation is then described by 



£(i) = max(0,2i?- \r x -f 2 \) 



(1) 



and the deformation rate £ (t) . Apart from material prop- 
erties, the dissipative and elastic components of the nor- 
mal force of particles in contact depend on the deforma- 
tion, the deformation rate, and the radii, 



F = (£, R X ,R 2 ) + F^ U, £, R 1 ,R 2 



(2) 



The functional form of these forces is model specific, see 
e.g. 0, [3, ||. Having specified the interaction forces, 
the dynamics of a ensemble of granular particles can be 
solved by a (force-based) Molecular Dynamics scheme. 

An alternative approach uses the concept of the coef- 
ficient of restitution, relating the normal component of a 
pair of particles before and after a collision, 



i'/i- 



(3) 



This concept does not consider the duration of a con- 
tact, that is, a collision is an instantaneous event. Con- 
sequently, it is assumed that the particles collide exclu- 
sively pairwise. This condition is justified if the mean 
flight time between collisions is much larger than the du- 
ration of a collision which restricts the range of applica- 
bility of the coefficient of restitution. The material prop- 
erties of the particles are, thus, assumed to assure short 
duration of contact and/or the particle number density 



of the system should be small enough (low collision fre- 
quency) to neglect multi-particle contacts. In practical 
applications, event-driven Molecular Dynamics simula- 
tions, based on the coefficient of restitution, deliver fre- 
quently satisfying results even for rather dense systems. 

Both concepts, interaction forces and the coefficient of 
restitution, can be applied to describe the dynamics of a 
granular system using either (force-based) Molecular Dy- 
namics or event-driven Molecular Dynamics. Describing 
the same physical systems, of course, the coefficient of 
restitution and the interaction forces must be closely re- 
lated. Indeed, integrating Newton's equation of motion 
for an isolated pair of particles colliding at time t = 0, 

mrt£ + F(£,t) =0; £(0)=0; £ = « (4) 

with m e ff = m\rnil (mi + m?) and 

v _ K(0)-^(0)]-[ri(0)-f 2 (0)] 

Rl + i?2 

to obtain the trajectory £(£), the coefficient of restitution 



(6) 



where t c is the duration of the collision. This computa- 
tion was performed for several interaction force models 
[E 0, EL HI- Albeit conceptually simple, even for simple 
force laws the algebra is rather technical. 

It is important that Eq. @ applies to particles in con- 
tact. Obviously, in the absence of adhesion, the interac- 
tion force between colliding granular particles is strictly 
repulsive. Formally, however, during the decompression 
phase where £ < the dissipative term F^ dls ^ in Eq. 
@ may overcompensate the pure repulsive conservative 
force F( el ) erroneously yielding an attractive total force, 
e.g. 0.0- 

In Molecular Dynamics simulations, therefore, the nor- 
mal force between particles is usually computed as F* = 
max(0, F), with F given in Eq. ([2]) which assures that 
only repulsive forces act. The force F* can, thus, be 
conveniently used in simulations. 
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The described artifact of negative interaction force 
originates from an inappropriate definition of the end of 
a collision at time f c . The duration of the collision, t c , 
however, is needed for the derivation of the coefficient of 
restitution by means of Eq. (J6]) . Whereas the beginning 
of a contact is well described by the condition £(0) = 0, 
the end of a collision at time t c is less trivial. 

For simplicity of the computation in the literature it 
was assumed that the end of a collision is determined by 
the condition 



particles we find 



f (t c ) = with t c > . 



(7) 



As described above, in the decompression phase it may 
happen that F K>£) < 0- This means the collision may 
be completed even before £ = 0. Thus, the surfaces of the 
particles lose contact slightly before the distance of their 
centers exceeds the sum of their radii. Consequently, 
the deformation of the particles may last longer than the 
time of contact and the particles gradually recover their 
spherical shape after they lost contact. The definition of 
the end of a collision 



F(t c ) = with t c > 



(8) 



takes the described scenario into account and assures that 
the particles interact exclusively repulsively. 

Obviously since erroneous attractive forces are ex- 
cluded by the improved condition for the end of collision, 
the resulting coefficient of restitution is expected to be 
larger for the definition Eq. © than the value obtained 
for the condition Eq. (J7J). 

Let us demonstrate the influence of the definition of t c 
to the coefficient of restitution for the simplest form of 
the interaction force, the linear dash-pot 



(9) 



Although neither the elastic nor the dissipative compo- 
nents are appropriate for the description of dissipatively 
colliding spheres (see below), the linear dash-pot model 
is frequently used in Molecular Dynamics simulations of 
granular systems. The main advantage of this model is 
the impact-velocity independent coefficient of restitution 
which follows from Eq. ((6|). Using the condition jjj), we 
obtain for the case of low damping (e.g. [HQ) 



exp 



/3tt 



and t c = 



(10) 



with u) = y/ujQ - (3 2 ; ujq = \/k/m c g; (3 = j/2m c g. Ob- 
viously, this result contradicts the assumption of non- 
attractive interaction since 



F 



(t c ) = F (t c ) , i (t c )) - F (0, -ev) = yev 



> 0. 



(11) 

For the condition Eq. ((5|) for t c , taking into account 
that there is only repulsive interaction between granular 



exp 



exp 



exp 
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/3e 



Ik."" 



(3 > w 



(12) 

It can be shown that the solutions, Eq. (|10[) and (fT2f are 
fundamentally different: for values of the parameter /3/ujq 
above one the duration of the collision t c diverges in case 
of Eq. (JTDJ) , that is, e = 0. Thus, the particles collide 
with finite velocity and stick together (dissipative cap- 
ture), despite our precondition of purely non-attractive 
interaction. The solution Eq. (fi"2")) does not reveal this 
unphysical behavior. For a detailed discussion see @]. 

The linear dash-pot model serves here only as an ex- 
ample to show that even for the simplest force laws the 
adequate characterization of the end of the collision mod- 
ifies the known results for the coefficient of restitution in 
a non-trivial way. For the case of the linear dash-pot, the 
definition of t c , Eqs. ([7]) or ©, changes the coefficient 
of restitution as a function of the material parameters k 
and 7 however, e is independent of the impact velocity v 
in both cases. 

It is the aim of this paper to compute the coefficient 
of restitution for the simplest physically consistent force 
law for viscoelastic spheres with regard to the definition 
Eq. ([5]) for t c . We will see that the appropriate choice 
of the condition for the end of the collision does not only 
change the dependence of the coefficient of restitution on 
the material parameters but also the functional form of 
its dependence on the impact velocity. 

As our main result we will show that for the definition 
of t c given by Eq. © the coefficient of restitution e 
is given by a series in powers of w 1 / 10 whereas for the 
definition of t c according to Eq. ([7]) e is a series in powers 
ofw 1 / 5 



II. VISCOELASTIC SPHERES 

We write the interaction force law for viscoelastic 
spheres Q as 



(13) 



The elastic part is given by the Hertz contact force [H 
with the elastic constant 



(14) 



3(1-^ 2 ) 



where Y is the Young modulus, v is the Poisson ra- 
tio and the effective radius of the colliding pair R c g = 
RiRj/ (Ri + Rj). The dissipative part, ~ V^S,, was de- 
rived independently in using different methods 
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but only the method in @ allows to derive the dissipative 
constant 



3 3i] 2 + 2r/i 



(1 - „ 2 ) (1 - 2^) 



Yi> 2 



(15) 



as a function of viscous material constants t\\i^ that re- 
late the dissipative stress tensor to the deformation rate 
tensor [l3| and the elastic constants Y and v. 

While the coefficient of restitution for the linear dash- 
pot model depends only on the material constants, it may 
be shown already from a dimension analysis that for vis- 
coelastic particles the coefficient of restitution cannot be 
independent of the impact velocity, [E EH EE EE E3 • It 
may be shown, moreover, either by scaling arguments [j| 
or in a more accurate way by a rather technical analy- 
sis [i| that the coefficient of restitution depends on the 
impact velocity as e = e (y 1 ^ 5 ) . The coefficient of resti- 
tution was obtained in for the definition © as a series 
expansion in powers of v 1 / 5 . (For an equivalent deriva- 
tion for viscoelastic discs see (18|.) In the following we 



derive the coefficient if restitution for the end of the col- 
lision given by Eq. ([5]). 



III. EQUATION OF MOTION 

Newton's equation of motion for the collision of vis- 
coelastic spheres reads 



with initial conditions 



'7 



£(0) = 0; £ = v. 



and the constants 

k = 



P 



_ 3 pA 
7 2 m cff 



(16) 



(17) 



(18) 



The natural unit of time is t s 



■ale 



k 2/5^ 1/5 which jg 



proportional to the duration of the undamped collision 
and the natural unit of length is £ sca ie = k~ 2 / 5 v 4 / 5 which 
is proportional to the maximal deformation. Adopting 
both natural units would reduce the number of free pa- 
rameters to one which reads 7fc~ 3 / 5 v 1 / 5 [5j. This indi- 
cates that the coefficient of restitution is a function of 
v 1 / 5 . For reasons which will become clear in the course 
of the computation (see explanation at Eq. (f3"Tj) ) it is 
not advisable to use the natural unit of length. Instead 
we adopt the length scale £ sca ic = fc~ 2 / 5 - We, thus, scale 
time and length as: 



t = 



£2/5 ' £2/5^1/5 

and arrive at the equation 

x + Qv^x^fx + v~ 2/5 x 3 / 2 = 
a;(0) = 0; i(0) = u 4/5 , 



(19) 



(20) 



where dots mean derivatives with respect to the scaled 
time t and 8 = 7 A; -3 / 5 . Note that the deformation £ or 
x are counted positive if the particles deform each other. 
The impact velocity £(0) or ±(0) has to be positive as its 
action increases the deformation. 



IV. TRAJECTORY 

First we have to determine the trajectory of the parti- 
cles during the collision. To this end we apply the method 
which was introduced in 0]. 

First we observe that the trajectory cannot be a series 
in integer powers of time due to the fact that the third 
and higher time derivatives of the deformation are sin- 
gular at x = 0. The deformation x — corresponds to 
the start of the collision and also to its end under the 
condition Eq. ([7]) . (Here we consider the collision for the 
condition Eq. ([8]), nevertheless, for the calculation we re- 
fer in several places to the end of the collision due to Eq. 
((7|) which we call the naive end of the collision.) As an 
example for such a divergence, the third time derivative 
of x reads: 



,(3) _ _. 



2v 2 / 5 



Xy X + 



2 
X - 



^3/5 



v 2/5- 



(3 



2V 1 / 5 ^fx 



(21) 



The last term diverges for x — > as for the beginning 
and the end of collision i ^ 0. It turns out that instead 
of integer powers the trajectory is a series of half-integer 
powers of r. The computation of the trajectory x(r) is 
explained in detail in appendix [X] The first few terms 
read 



x = v 
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(22) 



It turns out that this series converges very slowly which 
means that we need the series up to a high order (see be- 
low). The structure of this result becomes clear if we sort 
the terms in escalating powers of the damping parameter 
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P, The trajectory then takes the form 



X = | T 
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(23) 

The expressions in brackets do not contain any parame- 
ter except for pure numbers. They are, hence, universal 
functions which we shall call Xi(r), where the index i 
gives the power of P it is associated with. Note further- 
more that subsequent powers of r in each function differ 
by 5/2. The trajectory can be written compactly as 

x(t) = v 4/5 x (t) + pvxiir) + P 2 v 6/5 x 2 (t) + ... 



,4/5^ 
fe=0 



[pv 1 '^ x k {r) (24) 



The function xq is the trajectory of the undamped (P = 
0) collision. It is known [5J that it reaches its maximal 
compression at time 



I) 



21 [ — 

10 



1.609 . 



(25) 



The total duration of the undamped collision is t° = 
2 r max as the undamped trajectory is symmetrical with 
respect to the point of maximal compression. 

We proceed with computing the time of maximal com- 
pression of the damped problem along with the value of 
maximal compression. We use the Ansatz 



oo 

r max = r max + ^ t 
k=l 



a k p K v 



(26) 



and solve for the coefficients a k as explained in detail in 
Appendix [Bj The first coefficients a k are listed in Table 
HI The principal form of these and other similar expres- 
sions - power series in By 1 / 5 - can be derived by scaling 
arguments detailed in 5] . The maximal deformation can 
be obtained by Taylor expansion of Eq. 



(27) 



fe=0 



with the coefficients b k . We will not need these coeffi- 
cients explicitly, they can, nevertheless, be found in table 
III 



V. FINAL VELOCITY FOR THE NAIVE 
CONDITION 

Let us compute the final (naive) velocity, assuming 
the end of the collision according to Eq. ([7]). At first 
glance one might be tempted to compute the duration 
of collision with an Ansatz like r c = t® + 6t c and solve 
for the correction terms by performing a Taylor expan- 
sion around the undamped duration of collision. This 
method, however, fails due to the aforementioned singu- 
larity at x = 0. Instead we compute the final velocity 
indirectly: as we have an expression which is definitely 
valid for the first part up to the maximal compression 
we can construct the full solution by a kind of backward- 
shooting method. We start at the end of the collision 
where x = —v' (the final velocity v' being unknown yet) 
and let the time run backwards. The equation of motion 
for this inverse collision is identical to Eq. ([20]) 



(28) 



Pv'- 1/5 x invy / X — + v / - 2/5 x^ = 



£inv(0) = 0; x inv (0) 



except for the sign of the damping parameter P, since the 
inverse collision (in inverse time) is an accelerated colli- 
sion. Consequently, the trajectory of the inverse problem 
can be obtained from the solution of the direct collision, 



Eq. (|22j) . by simply substituting P — > — p and v 
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+ O (r 6 ) (29) 



The same is true for the maximal compression of the 
inverse collision, 



xZ x = v'*/ 5 '£(-l) k b k p k v' k / 5 , (30) 



with the same numerical coefficients b k as in Eq. ([27)1 . 

As the inverse collision problem is just a reformulation 
for the original collision problem both maximal compres- 
sions have to be the same, 

%max — ^max i (31) 

which is an equation for v'. From these arguments the 
choice of our length scale, Eq. (fl~9|) . becomes evident: 
choosing the natural unit of length, fc -2 / 5 ?; -1 / 5 , the di- 
rect and the inverse collision problem would have dif- 
ferent length scales as the initial velocity of the inverse 
collision is v' ^ v. 

In order to solve Eq. (|3"Tj) for v' we use the Ansatz 



v + PSv! + P 2 6v 2 + 



(32) 
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and solve for the corrections Svi . Using the definition Eq. 
© this yields the coefficient of restitution of the form 



e(vY 



1 + dpv 1 ' 5 + c 2 /?V/ 5 + 



(33) 



Note that we determined the final velocity v' at x = 0, 
that is, this result for e(v) corresponds to the condition 
Eq. ([7]) for the end of the collision. Based on the trajec- 
tory derived so far, in the next section we will derive the 
coefficient of restitution that corresponds to Eq. ([8]) . 

The calculation of the coefficients Cfc in Eq. (f3"3"|) is 
explained in Appendix [B] the numerical values of the 
first coefficients are shown in Table HI 



VI. PREMATURE END OF THE COLLISION 

Up to here we calculated the solution of the equation 
of motion, Eq. fT6]) . in the interval (£ = 0, £ = v) (start 
of the collision) to (£ = 0, £ = v') (end of the collision) 
or the scaled Equation (|2"0)) in the corresponding interval 
x = in the beginning and x — in the end, respectively. 
The velocity at the end of this trajectory, v' , led us to the 
coefficient of restitution corresponding to the condition 
Eq. ©. 

As discussed before, however, the velocity v' corre- 
sponds to a negative interaction force, in contradiction 
to the purely repulsive interaction of viscoelastic gran- 
ular particles. Therefore, the collision does not end at 
x = but before, when the interaction force becomes 
zero. This condition corresponds to the condition Eq. 
©■ 

We take this premature end of collision into account 
and, thus, look for the earliest point in time T during the 
inverse collision when the acceleration vanishes. Setting 
i inv = in Eq. (J^SJ) yields 



(3v n / 5 x inv {T) =x inv (T). 



(34) 



For small (3v' 1 / 5 we obtain T to lowest order by approx- 
imating Xinv by i/ 4 / 5 r which yields 



T^f3v' 1/5 . 
The solution to higher order reads: 



T = /V 1 / 5 + —pV 2 v'V w + A/?v 6 / 5 
' 35 H 75 M 



(35) 



21271 
2734875 



^17/2^17/10 



(36) 



The details of this calculation can be reviewed in Ap- 
pendix [Cj The value of x- lnv at this point in time is 



1 + tV ,5/v1,2 + ^ v 



Going back to the original units of time we obtain the 
final velocity for the case of the condition Eq. ([H]), 



final 



1 + ±06/2 v n/2 + 

15 H 21CT 



(38) 



Inserting the expression for v' one arrives at the final 
solution 

e = l-1.153/% 1/5 + 0.798/3V /5 + 0.267/3 5/ V /2 + ... 

OC 

= l + J2hk(3 k/2 v k / w (39) 

The details of this computation are shown in Appendix 
[Cl The coefficients hk are pure numbers; the first 20 of 
them can be found (to a higher precision as in the expres- 
sion above) in table [Til As the coefficient of restitution 
e only depends on f3v 1 ^ 5 (including half powers of this 
term) we show the velocity dependence in this universal 
form in Fig. [TJ 

The analytical results, Eqs. (|33p and (|3"9")1 . are com- 
pared with the numerical solution of the equation of mo- 
tion (|20p . In the interval shown in Fig. [TJ the analyti- 
cal results agree with the numerical results almost per- 
fectly. Beyond the shown interval the solutions start to 
deviate. As an example in physical units we consider a 
sphere that collides with e = 0.8 at v = lm/sec, e.g. 
a rubber sphere. By numerically solving Eq. (|39|) we 
obtain j3 = 0.2 sec 1 ^ 5 /m 1 / 5 . Consequently, the range of 
velocity shown in Fig. [U j3v x / b < 0.3 corresponds to 
v < 7.5m/sec. From the good agreement between the 
analytical and numerical solutions in this interval we con- 
clude that the range of validity of the solution, Eq. ([39)l . 
is at least v < 7.5m/sec. For materials with smaller 
damping constant /3 the range of validity is larger. 

Albeit in Fig. [ljnumerical and analytical results almost 
coincide we note that the deviation for the improved con- 
dition, Eq. ([5]), exceeds the deviation for the naive con- 
dition by several orders of magnitude. This can be seen 
from the coefficients hk (see table HTJ which decrease only 
slowly for increasing k. Thus, to obtain a good precision 
for fiv 1 / 5 close to unity requires a very large number of 
coefficients hk- 

For large velocities or large damping both velocity de- 
pendencies, Eqs. ([33| and (|39p . reveal a remarkable dif- 
ference: For the naive condition, Eq. ([7]), the coefficient 



of restitution decays asymptotically as e 



For the 



improved condition, Eq. ([S]), the asymptotics is compat- 



ible with a power law of e 
are shown in Fig. [21 



-0.331 



Both asymptotics 



(37) 



VII. CONCLUSION 

We described the collision of a pair of particles which 
interact repulsively according to the force law, Eq. fjjjj, 
valid for viscoelastic spheres. In a physically consistent 
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FIG. 1: The velocity dependence e(v) for both conditions for 
the end of the collision, Eq. (0 (naive) and Eq. Q (im- 
proved). The upper panel shows the dependence on velocity, 
the lower panel shows the dependence on (iv x ^ . For both 
panels the interval shown in the abscissa is (almost) equiva- 
lent. The numerical solution of Newton's equation of motion, 
Eq. (|20[) . agrees almost perfectly with the analytical curves 
shown here. Beyond the shown interval there are increasing 
discrepancies between theory and simulation. 



description the end of the collision is determined by the 
instant during the expansion when the interaction force 
vanishes, = 0, (a) but not by the naive condition 
=0 (b) which corresponds to the instant when the 
distance of the centers of the particles coincides with the 
sum of their radii. This becomes obvious when looking at 
the interaction force at the end of the collision: For con- 
dition (b) the interaction force becomes attractive which 
contradicts the precondition of purely repulsive interac- 
tion. The reason for this behavior is the delayed recovery 
of the particles, that is, the surfaces of the particles lose 
contact already slightly before the compressed particles 
recovered their spherical shape. 

The choice of the condition for the end of the collision, 
(a) or (b), has a drastic effect on the resulting velocity 
dependence e(v) of the coefficient of normal restitution. 
Instead of a series in u 1 / 5 obtained for the naive condition 



FIG. 2: Asymptotics for e(v) for the naive end-of-collision 
condition, Eq. 0, (top) and the improved condition, Eq. 
(f8jl. (bottom). Both are representable by a simple power law. 
For the naive condition the exponent is —1, for the correct 
condition it is close to —1/3. 

(b) 0, [f|, for the physically consistent end-of-collision 
condition (a) we obtain a series in u 1 / 10 where the odd 
powers of w 1 / 10 are solely due to the end-of-collision rule. 
The analytical results agree almost perfectly with the 
numerical integration of Newton's equation of motion for 
colliding viscoelastic spheres. 

We evaluated the result for e(v) for realistic material 
properties for the cases (a) and (b) and obtained a no- 
ticeable difference of up to about 20%, depending on the 
material properties. The range of validity of our result 
was estimated by about 10 m/sec for a soft, rather dis- 
sipative material such as rubber. For a more elastic ma- 
terial, corresponding to a larger coefficient of restitution, 
the range of validity is significantly larger. Our analytical 
results deviate from the numerical results for /Ju 1 / 5 > 0.9 
which may be attributed to the properties of the series, 
Eq. (|39| . which converges slowly for large (3v 1 / b and 
whose convergence is not even clear for [iv 1 ^ > 1. 

For large impact velocity we can, however, still obtain 
numerical results which reveal another drastic difference 
between the conditions (a) and (b): For both conditions, 
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asymptotically e(v) follows a power law. For the naive 



condition (b), however, we obtain e 



whereas for 



the physically consistent condition (a) we find e ~ v~ x / 3 . 

The influence of the end-of-collision condition on the 
coefficient of restitution for viscoelastic particles is in 
marked contrast to the corresponding result obtained for 
the linear dash-pot model [g]. Here the choice of the 
condition (a) or (b) would result a modified coefficient 
of restitution which is, nevertheless, independent of the 
impact velocity in both cases. 
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APPENDIX A: COMPUTATION OF THE 
TRAJECTORY 

Equation ([20| for the trajectory x(t) of the particles' 
relative motion in the scaled variables x and r, 



x + (3v~ 1/5 xy^x + v~ 2/5 x 3/2 = 



x(0) = 0; x{0) = v 



(Al) 



is solved by series expansion. As explained in the text, 
an expansion in powers of t fails, instead we expand in 
powers of y^r. Using the Ansatz 



x(t) = w 4/5 t [1 + R{t)} 



(A2) 



Eq. (|Aip turns into 

2R + tR + f3v 1/5 r 1/2 (l + R + tR) VTTR 
+ r 3/2 [l + i?] 3/2 = 



(A3) 



R(0)=Q; i?(0) = 0. 

The term R may be expanded in powers of i/t, 

R(t) =a + cut 1 / 2 + a 2 T X + a 3 r 3 / 2 . . . (A4) 

Inserting Eq. (|A4[) into Eq. (|A1[) and comparing equal 
powers of t 1 / 2 we find a$ = ai = a 2 = 0, that is, the first 
non-trivial contribution is C(r 3 / 2 ). This fact simplifies 
the subsequent computer algebra considerably. 



We determine the coefficients a 3 , 04, . . . in escalat- 
ing order using an iterative procedure. In the first step 
we determine a 3 while <Zj (i > 3) stay undetermined. 
The corresponding term for R of the order 3 is denomi- 
nated by i?3 = a 3 r 3 / 2 + OIt 2 ) the next order is 4 with 
i?4 = a4T 2 + 0(t 5 / 2 ), etc. In other words, Ri contains 
all contributions of order 0(W 2 ) and higher. In each 
step i of the iteration we derive a differential equation 
Gi[Ri, Ri,Ri] — for Ri. 

We demonstrate the procedure for the first terms of a 
series up to the term agT 9 / 2 . For the first step, i = 3, 
we expand (1 + R) 1 / 2 and (1 + Rf/ 2 in Eq. flM} up to 
the necessary order for R3. Since N = 9 and the lowest 
order of t in R is 3, we need the expansion up to the 
third term, 



\/l + R 3 =1 



2 



Ri 



P3 

■"-3 

16 



/, „ \S/9 1 3i? 3 3i?? Rq 

(1 + R3Y' 2 =1 H H - 

y iJ 2 8 16 



(A5) 



Equation (|A3j) reads then 

G 3 [R 3 ] = 2R 3 +tR 3 + /3v 1 ^t 1 ^ 2 x 

x(l + ^R 3 + -#3 + tR 3 + ^tR 3 R 3 

+ - 3/2 (l + ^3 + ^)=0, (A6) 



where terms of order r 10 / 2 and higher are neglected. 

The desired coefficient 03 is now isolated by the formal 
transformation 



R 3 = a 3 r 3 / 2 + Ri 



(A7) 



which establishes the first iteration step. In general, we 
replace 



Ri — a-iT 1 / 2 + Ri+i 



(A8) 



insert this into Gi(Ri) — where only terms of rele- 
vant order are taken into account. Then we consider the 
term O (W 2-1 ) and determine a*. After substituting 
back into Gi we are left with the next order equation 
Gi+i (i?i+i) = which is then solved in the same way, 
etc. 

We insert Eq. (|A7f into Eq. and obtain 



15 

2i? 4 + tR a + ( — a 3 + f3v 



1/5 ) t 1 ' 2 + r 3 / 2 + 3Pv^a 3 T 2 + ^a 3 T 3 + l(3v^a 2 3 T 7 ' 2 + |a 2 r 9 / 2 
J 2 8 8 



hv^RiT 1 ' 2 + Irat 3 ' 2 + /V/ 5 i? 4 r 3/2 + jW /5 a 3 i? 4 T 2 + l^^Rlr 1 ' 2 = (A9) 
2 2 2 8 



8 

where again terms of irrelevant order were skipped. The terms in brackets of lowest order 1/2 (i/2 — 1 in general) 
allows for the computation of the first non-trivial coefficient 03 = — (4/15)/9t; 1 / 5 . We insert 03 into Eq. (|A9p to obtain 
the next equation for the computation of 04: 

G 4 [Rt] = 2R, + tR 4 + Ipv^R^ 2 + Uv^Rlr 1 ' 2 + r 3 ' 2 + ^t 3 ' 2 + f3v^R^ 2 - l/3 2 v 2 / 5 R 4 r 2 

_ 1^ W 2/5 T 2 _ 2^1/5 T 3 + A /3 3 u 3/5 T 7/2 + J. ^2/5^/2 = q (AlQ) 

5 5 25 75 

The next iteration step R 4 — a 4 r 2 + R$ leads to 

2i? 5 + tR 5 + 6a 4 r + t 3 ' 2 - U 2 v 2 I*t 2 + 1/^/5^5/2 _ 2^i/5 T 3 + (3 + 2 ^ 8/5 \ t?/2 

5 2 5 \2 25 / 

+ hv^R^ 2 + ^ 5T 3/2 + /5 w l/5^ 5r 3/2 _ ^2/5^2 _ ^2/5^4 + ^1/5^9/2 + 2 ^2/5^/2 _ Q 

2 2 5 5 8 T5 

(All) 

From the terms of lowest order we find = 6a 4 r, that is a 4 = 0. We insert this into Eq. (I Al 1|) : 

G 5 [R 5 ] = 2R 5 + tR 5 + r 3 / 2 - -f} 2 v 2 l\ 2 - ^v^t 3 + A/3V/ 5 t 7 / 2 + hy^R^ 1 ' 2 + Ir 5 t 3 ' 2 

5 5 25 2 2 

+ /^/^sr 3 / 2 - -/3 2 i; 2 / 5 i? 5 r 2 + A/3V/V 9 / 2 = . (A12) 
5 75 

With i?5 = (Z5T 5 / 2 + i?6 the latter equation turns into 

2R 6 + tR 6 + (^a 5 + l) r 3 / 2 - -^ 2 v 2 '\ 2 + _ ^i/s^j 

+ ^ 5 r 4 + (-l/3 2 v 2 ^a 5 + ^f3 2 v 2 A r 9 / 2 + hv^R^r 1 ' 2 + /V/ 5 i? 6 r 3 / 2 + |i? 6 r 3 / 2 = . (A13) 
2 \ 5 75 / 2 2 

From the lowest order terms O (t 3 / 2 ) we obtain 05 = —4/35. We insert 

G 6 [R e ] ee 2R e + tR 6 - U 2 v 2 /\ 2 - Uv^r 3 + ^v^r 7 ' 2 - ^ + ^V/V 9 / 2 + hv^R^ 2 
5 7 25 35 525 2 

+ *-R 6 t 3 ' 2 + pv^R^' 2 = (A14) 

iterate, i?6 — o,qt 3 + i?7, and obtain 

2R 7 + R7T + (l2a 6 ^ 2 v 2 ^ r 2 - ^v^r 3 + Q/V^ + |/?V/ 5 ) r 7 / 2 - |r 4 



T 3 + 2 /3 3 u 3/5 T 7/2 



{r e + Jl^ 275 ) r9/2 + ^ 1/5i?7rl/2 = • ( A15 ) 



From the term of lowest order we find o,q = j3 2 v 2 / 5 /15. We insert a§ into Eq. (|A15p for the next order equation 

G 7 [R 7 ] ee 2R 7 + rR 7 Uv^r 3 + ^-B 3 v 3 '\ 7 l 2 - £-r* + ^P 2 v 2/5 r 9 / 2 + h v ^R 7 ^ 2 = . (A16) 
7 50 35 1050 2 

Iterating R 7 = a 7 r' 7 1 2 + R$ yields 

2R S + R s r + ^a 7 r 5 / 2 - + 12^3/5^/2 + 3^,1/5 4 _ 6 ^ + ^1^2/^9/2 + = Q (A17) 

4 7 50 2 35 1050 2 

and from lowest order ~ t 5 / 2 we obtain: a 7 = 0. We insert this solution, and substitute Rs = asT 4 + i?g: 

2R9 + tR 9 + (20a 8 - Q -(3v^ r 3 + ^WS' 2 - |r 4 + ^v^a 8 + ^/3 2 - 2/5 ) r»' 2 = (A18) 

I 



We insert the solution as = j^flv 1 ^: replace Rg = agr 9 / 2 + R 



in. 



G 9 [R 9 ] = 2R 9 + rR 9 + ^/3V/ V/ 2 - Ar 4 2i? 10 + R 10 r + (^a 9 + ^P 3 v 3/ A r 7 / 2 

50 35 \ 4 50 / 

+ ^ V/V9/2 = ( M9 ) _l T 4 + A^2/5 r 9/2 = (A20) 

300 35 300 M v ; 
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and obtain ag = — (38/2475)/3 3 u 3 / 2 . Inserting this solu- 
tion and substituting R w = aiQT° + i?n yields 

2R U + r Rll + ( 3 0a 10 ~ |) r 4 + ^/3V /5 t 9 / 2 = 

(A21) 

and, thus, aio = 1/175 which is the last coefficient which 
can be obtained from the expansion, Eq. (| A5|l 

To achieve an acceptable accuracy of the final result, 
the expansion Eq. has to be performed up to high 
orders in fiv 1 / 5 . To accurately compute the necessary co- 
efficients hk one needs accurate functions Xk of the same 
index k. For the chosen accuracy (20 coefficients hk) the 
expansion of the trajectory has to be performed up to 
an order as large as 150. We employ computer algebra 
(maple) which turns the described algorithm into only 
a few lines of code. For the computation we abbreviate 
A = /3d 1 / 5 , s = y/r, Rd and Rdd stand for dR/dr and 
d 2 R/dr 2 , and N is the order of the expansion. 

restart ; 
N := 150; 

dgl : =2*Rd+s~2*Rdd+(A*s+s~3) * (1+R) * (3/2) 

+A*s~3*Rd*sqrt(l+R) ; 
dgl:=convert(taylor(dgl,R,N) ,polynom) : 
solution :=0; 
for i from 3 to N do 

dgl : =subs (Rdd= (i* (i-2) /4) *a*s~ (i-4) +Rdd , 

Rd=(i/2)*a*s~ (i-2)+Rd,R=a*s~i+R,dgl) : 
dgl:=mtaylor(dgl, [R,s,Rd,Rdd] ,N, [i,l,i,i]) : 
tmp : =expand(coef f (coef f (coef f (dgl ,R,0) , 

Rdd.O) ,Rd,0)) ; 
asol:=solve(coeff (tmp, s, i-2) ,a) : 
print (i , asol) : 

dgl : =simplif y (subs (a=asol ,dgl) ) : 
solution : =solution+asol*s~i : 
end do : 

solution:=v~ (4/5) *s~2* (1+solution) : 
solution: =subs (A=beta*v~ (1/5) , solution) : 
f out :=f open (" ./solution" , WRITE) ; 
f printf (f out , "°/a\n" , solution) ; 
f close (f out) ; 

APPENDIX B: TIME AND VALUE OF 
MAXIMAL COMPRESSION AND THE SERIES 

C naive 

( ,! ) 

The first ingredient for the actual computation of 
£naive(^) is the maximum compression. To this end, we 
first compute at which time r max this maximum compres- 
sion is achieved. The time of maximum compression will 
be determined by Taylor-expansion of the expression 

iftL* + H=°- ( B1 ) 

Here the time T" lax of maximum compression of the un- 
damped collision as given by Eq. (I25p is taken as a refer- 
ence. In terms of the universal functions introduced 



in Eq. (|23|) the Taylor expansion takes the form 

nf3 n f> r lk+l r k 

i (tL, + Sr) = E ( B2 ) 

which motivates the representation of St as a series of 
the form 

rip 

St = Y^ a n P n v n/5 . (B3) 

n=l 

We insert Eq. (|B3|) into the Taylor expansion, collect 
coefficients in powers of (3 and solve successively for a n . 
The result is shown in Tab. HI 

In the same way the maximum compression a; max 
can be computed by performing the Taylor-expansion of 
x ( T max + which is of the form 

^ d k T ST k 

i=0 k=0 

suggesting the series 

rib 

x max = v 4 / 5 J2 b n (3 n v n/5 (B5) 

n=0 

The coefficients b n are shown as well in Tab. [T] The first 
coefficient bo is the maximum compression for the un- 
damped problem. To actually compute the coefficient of 
normal restitution without regard of the premature loss 
of contact we have to match the maximum compression 
of the direct and the inverse collision, i.e. we have to 
solve 

= <L («') (B6) 

for v' with 

rip 

x max (v) = v^J2 b ^ nvn/5 ( B7 ) 

n=0 
rig 

x max (v') = x/ 4 / 5 ]T(-ir& n /?V"/ 5 (B8) 
n=0 

In the Maple program the function x max (v) is called h(v) , 
the function x max (v) is called hm(v). Using the Ansatz 

v' = v(l + CnP n v n/5 ^J = ve nilivc (v) (B9) 

we can solve for c„ by expanding the expression Eq. (|B6[) 
for small (3 and collect orders. The first Ck are shown in 
Tab. H 

restart; Digits := 20; 
nb := 20; 

fin := fopen(" ./solution" , READ): 
xin := fscanf (fin, "7.a") : 
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TABLE I: The first numerical coefficients of the expansions 
Eq. dB3j 







bi 




o 




1 093362074 




1 


-0.2867471220 


-0.5044548926 


-1.153448854 


2 


0.1048589922 


0.2840430192 


0.7982665553 


3 


-0.04868406400 


-0.1702207776 


-0.5228825609 


4 


0.02543116890 


0.1055007088 


0.3487426678 


5 


-0.01423658282 


-0.06684871371 


-0.2330981260 


6 


0.008337660013 


0.04303949229 


0.1566821477 


7 


-0.005039737366 


-0.02805108430 


-0.1058187828 


8 


0.003118137108 


0.01846085121 


0.07176528242 


9 


-0.001964027745 


-0.01224618562 


-0.04885717237 


10 


0.001254701962 


0.008177589114 


0.03337347194 



x := simplify (subs (s = sqrt(t), xin[l])): 
tchalf := (4/5)~(3/5)*GAMMA(2/5)*GAMMA(l/2)/ 

(2*GAMMA(9/10)) : 
x := subs(t = tchalf +dt, x) : 
xdot := evalf (taylor (dif f (x, dt) ,dt=0,nb)) : 
xdot := convert (xdot , polynom) : 
dt := sum(a[ , i']*beta~'i , *v~((l/5)* , i'), 

' i ' = 1 . . nb) : 
for i to nb do 

a[i]:= solve (coeff (xdot , beta, i) , a[i]) 
od: 

hh := convert (evalf (taylor (x, beta, nb+D), 

polynom) : 
xmax := unapply(hh, v) : 

xmaxinv := unapply (subs (beta = -beta, hh) , v) : 
u := v*(l+sum(c['k']*beta~'k'*v~((l/5)*'k'), 

'k' = 1 . . nb)) : 
d := convert (taylor (xmax (v) -xmaxinv (u) , 

beta, nb+1) , polynom) : 

for i to nb do 

c[i]:= solve (coeff (d, beta, i),c[i]) 
od: 

fout := fopen(" ./coefficients" , WRITE): 
for i to nb do 

fprintf (fout, "°/„a\n", c [i] ) 
od: 

f close (fout) : 



APPENDIX C: PREMATURE LOSS OF 
CONTACT 



As the moment of actual loss of contact is close to the 
naive end of contact we will use the inverse collision to 
compute the time T and velocity v" at loss of contact. 
Using the condition £j nv = we obtain the equation for 
T: 



Approximating x; nv as u' 4 / 5 T we obtain the leading order 
of 



(C2) 



T = Pv n / 5 



After canceling the common prefactor v' 4 / 5 Eq. (|C1|) 
does only depend on the combination dv' 1 ' 5 . Therefore, 
one can easily guess the principal form of T: 



fc/5 



(C3) 



k=2 



Inserting this Ansatz into Eq. (|C1[) , collecting orders and 
solving for dk yields 



T = Bv' 1 / 5 + — /? 7 /V 7 / 10 + — /3V 6 / 5 
35 75 



21271 
+ 2734875 

The final solution now reads 



017/2^/17/10 



(C4) 



V"=v(l + + —(3 5 V + J^L /3 15/2 w 3/2 

V 15 M 150 M 160875^ 



453263 



3 w v 2 + 



40540500' 

Inserting the known solution for v' we obtain 



v" = v(l + J2h k P k/2 v k / w 

\ fc=0 / 
oo 

e(v) = l + J2hkP k 



1 k/2 y k/10 



(C5) 

(C6) 
(C7) 



fc=0 



8v n / 5 x inv (T) =x inv (T) 



(CI) 



The first values of are tabulated in Table HT1 

restart : Order : =20 : nd : =4 : Digits : =20 : 
fin:=fopen(" . / solution" , READ) : 
L : =f scanf (fin, ""/a") : f close (fin) : 
x:=convert (taylor (L[l] , s, Order) , polynom) : 
xinv:=subs(s=sqrt(T) ,subs(beta=-B*B,x)) : 
xinvdot : =dif f (xinv,T) : 

eqn: =simplif y (xinv-B*B*v~ (l/5)*xinvdot) : 
T:=B"2*v"(l/5)*sum(d['k']*B~(5*'k , )*v"('k , /2) , 
'k'=0. .nd) : 

eqn: =expand(eqn) : 
eqn:=series(eqn,B,2*0rder+l) : 
for i from to nd do 

d[i] :=solve(coeff (eqn,B,2+5*i) ,d[i] ) : 
od: 

vpp : =convert (series (v~ (1/5) *xinvdot ,B, 2*0rder+l) , 

polynom) : 
f in:=f open(" . /coefficients" ,READ) : 
for i from 1 to Order do 

L : =f scanf (fin, "°/a") : 

c[i] :=L[1] : 
od: 

f close (fin) : 
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TABLE II: The first numerical coefficients of the expansion 
(|C6|) . The coefficients /12 and /14 are identical to the first 
coefficients in the original expansion e(«), i.e. Y12 = ci and 
h 4 = c 2 . 



£ 


hi 





1 


1 





2 


-1.153448856 


3 





4 


0.7982665581 


5 


0.2666666667 


6 


-0.5228825657 


7 


-0.4613795424 


8 


0.3487426737 


9 


0.4523510496 


10 


-0.1464314644 


11 


-0.3677282992 


12 


-0.0432489833 


13 


0.2818042325 


14 


0.1478525872 


15 


-0.1794420590 


16 


-0.1784660326 


17 


0.06593358882 


18 


0.1713586178 


19 


0.0252498223 


20 


-0.1379234986 



vprime:=v*(l+sum(c['k']*B~(2*'k')*v~('kV5) , 

'k'=l. .Order)) : 
vpp : =convert (series (subs (v=vprime , vpp) , 
B,2*0rder+1) , 
polynom) : 
epsilon:=simplify(vpp/v) ; 
f out :=f open ("hk", WRITE) : 
for i from 1 to 2*0rder do 

h[i] : =simplif y (coef f (epsilon.B , i) /v" (i/10) ) : 
fprintf (fout,"%a,\n",h[i]) : 
od: 
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